Dicom loading and visualization¶
Cargar Archivos DICOM¶
Cargo los archivos DICOM desde el directorio especificado. Trabajo inicialmente con un subconjunto de los archivos para testear.
¿Con que CT trabajo? No me ha quedado claro qué subconjunto de archivos usar!
- ct_directory = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/03-14-1998-NA-CT ABDOMEN WWO CONT-48924/4.000000-Recon 2 3 PHASE LIVER ABD-87008"
- ct_directory_2 = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/03-14-1998-NA-CT ABDOMEN WWO CONT-48924/2.000000-PRE LIVER-48636"
Libraries¶
In [ ]:
import os
import pydicom
import numpy as np
import matplotlib.pyplot as plt
from typing import List
import cv2
import matplotlib.colors
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.animation import FuncAnimation
Aqui alojamos las funciones¶
In [ ]:
def load_dicom_files(directory: str, file_extension: str = '.dcm'):
"""Load DICOM files from a specified directory that match the given file extension."""
return [pydicom.dcmread(os.path.join(directory, f)) for f in os.listdir(directory) if f.endswith(file_extension)]
def sort_dicom_files_by_instance_number(dicom_files: List[pydicom.Dataset]) -> List[pydicom.Dataset]:
"""Sort DICOM files by the Instance Number."""
return sorted(dicom_files, key=lambda x: int(x.InstanceNumber))
def sort_dicom_files_by_position(dicom_files):
"""Ordena archivos DICOM según la posición del paciente en la imagen."""
return sorted(dicom_files, key=lambda x: float(x.ImagePositionPatient[2]))
def sort_dicom_files(dicom_files):
"""Sort DICOM files by available headers. Handles both CT and segmentation files."""
def get_sort_key(file):
# Define sorting keys: AcquisitionNumber if available, otherwise ImagePositionPatient
acquisition_number = getattr(file, 'AcquisitionNumber', None)
position_patient = getattr(file, 'ImagePositionPatient', None)
# Convert ImagePositionPatient to a sortable form if available
if position_patient:
position_patient = float(position_patient[2])
return (int(acquisition_number) if acquisition_number else 0, position_patient if position_patient else 0)
'''
Primero ordena por AcquisitionNumber y, dentro de cada número de adquisición,
ordena por la posición del paciente a lo largo del eje Z (ImagePositionPatient[2]).'''
return sorted(dicom_files, key=get_sort_key)
def check_acquisition_uniformity(dicom_files: List[pydicom.Dataset]):
"""Check if all DICOM files have the same 'Acquisition Number'."""
acquisition_numbers = {img.AcquisitionNumber for img in dicom_files}
if len(acquisition_numbers) > 1:
print("Multiple acquisition numbers found. Ensure all CT images are from a single acquisition.")
def get_pixel_arrays(dicom_files: List[pydicom.Dataset]) -> np.ndarray:
"""Extract pixel arrays from sorted DICOM files and stack them into a single numpy array."""
return np.stack([file.pixel_array for file in dicom_files])
def display_images(ct_image, seg_image):
"""Display CT and Segmentation images side-by-side."""
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
plt.imshow(ct_image, cmap='gray', aspect='auto')
plt.title('CT Slice')
plt.subplot(1, 2, 2)
plt.imshow(seg_image, cmap='gray', aspect='auto')
plt.title('Segmentation Slice')
plt.show()
def create_mip(ct_volume: np.ndarray):
"""Create and visualize Maximum Intensity Projection (MIP) of a CT volume."""
mip_axial = np.max(ct_volume, axis=0)
mip_coronal = np.max(ct_volume, axis=1)
mip_sagittal = np.max(ct_volume, axis=2)
projections = [mip_axial, mip_coronal, mip_sagittal]
titles = ['MIP Axial', 'MIP Coronal', 'MIP Sagittal']
plt.figure(figsize=(15, 5))
for i, mip in enumerate(projections):
plt.subplot(1, 3, i + 1)
vmin, vmax = np.percentile(mip, [2, 98])
plt.imshow(mip, cmap='gray', aspect='auto', vmin=vmin, vmax=vmax)
plt.title(titles[i])
plt.axis('off')
plt.show()
def visualize_planes(ct_volume: np.ndarray):
"""Visualize axial, coronal, and sagittal planes of a CT volume."""
coronal_plane = ct_volume[:, :, ct_volume.shape[2] // 2]
sagittal_plane = ct_volume[:, ct_volume.shape[1] // 2, :]
axial_plane = ct_volume[ct_volume.shape[0] // 2, :, :]
plt.figure(figsize=(15, 5))
for i, (plane, title) in enumerate(zip([axial_plane, coronal_plane, sagittal_plane],
['Axial Plane', 'Coronal Plane', 'Sagittal Plane'])):
plt.subplot(1, 3, i+1)
plt.imshow(plane, cmap='gray', aspect='auto')
plt.title(title)
plt.axis('off')
plt.show()
def apply_window(image: np.ndarray, window_center: int, window_width: int) -> np.ndarray:
"""Apply window leveling to an image based on the given center and width."""
lower_bound = window_center - 0.5 * window_width
upper_bound = window_center + 0.5 * window_width
return np.clip((image - lower_bound) / window_width, 0, 1)
In [ ]:
# Data paths
ct_directory = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/03-14-1998-NA-CT ABDOMEN WWO CONT-48924/4.000000-Recon 2 3 PHASE LIVER ABD-87008"
seg_directory = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/12-27-1997-NA-AP LIVER PRT WWO-67834/300.000000-Segmentation-39839/1-1.dcm"
#ct_directory_2 = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/HCC_007/manifest-1643035385102/HCC-TACE-Seg/HCC_007/03-14-1998-NA-CT ABDOMEN WWO CONT-48924/2.000000-PRE LIVER-48636"
Cargamos las imagenes y las ordenamos¶
Cargamos y ordenamos las imagenes CT
In [ ]:
# Load and process DICOM CT images
ct_files = load_dicom_files(ct_directory)
sorted_ct_files = sort_dicom_files(ct_files)
Cargamos la segmentation image
In [ ]:
# Load and process segmentation image
segmentation_image = pydicom.dcmread(seg_directory).pixel_array
sorted_seg_img = sort_dicom_files(segmentation_image)
mid_slice_index = segmentation_image.shape[0] // 2
Simple visualización para saber que todo ha ido bien
In [ ]:
# Visualize a slice of CT and Segmentation images
ct_images = get_pixel_arrays(sorted_ct_files)
display_images(ct_images[8], segmentation_image[8])
Estudio de encabezados¶
In [ ]:
# Print details and dimensions
print("Details of CT images sorted by Image Position Patient:")
for img in sorted_ct_files:
print(f"Acquisition Number: {img.AcquisitionNumber}, Image Position Patient: {img.ImagePositionPatient}")
print("Dimensiones del arreglo de segmentación:", segmentation_image.shape)
print("Dimensiones del arreglo de imágenes CT:", ct_images.shape)
Details of CT images sorted by Image Position Patient: Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -260.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -257.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -255.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -252.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -250.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -247.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -245.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -242.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -240.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -237.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -235.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -232.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -230.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -227.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -225.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -222.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -220.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -217.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -215.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -212.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -210.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -207.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -205.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -202.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -200.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -197.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -195.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -192.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -190.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -187.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -185.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -182.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -180.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -177.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -175.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -172.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -170.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -167.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -165.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -162.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -160.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -157.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -155.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -152.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -150.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -147.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -145.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -142.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -140.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -137.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -135.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -132.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -130.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -127.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -125.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -122.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -120.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -117.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -115.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -112.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -110.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -107.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -105.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -102.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -100.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -97.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -95.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -92.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -90.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -87.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -85.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -82.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -80.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -77.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -75.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -72.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -70.000000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -67.500000] Acquisition Number: 1, Image Position Patient: [-199.800003, -180.000000, -65.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -350.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -347.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -345.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -342.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -340.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -337.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -335.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -332.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -330.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -327.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -325.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -322.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -320.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -317.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -315.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -312.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -310.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -307.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -305.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -302.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -300.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -297.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -295.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -292.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -290.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -287.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -285.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -282.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -280.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -277.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -275.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -272.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -270.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -267.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -265.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -262.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -260.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -257.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -255.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -252.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -250.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -247.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -245.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -242.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -240.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -237.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -235.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -232.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -230.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -227.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -225.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -222.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -220.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -217.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -215.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -212.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -210.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -207.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -205.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -202.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -200.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -197.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -195.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -192.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -190.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -187.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -185.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -182.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -180.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -177.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -175.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -172.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -170.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -167.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -165.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -162.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -160.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -157.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -155.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -152.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -150.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -147.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -145.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -142.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -140.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -137.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -135.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -132.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -130.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -127.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -125.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -122.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -120.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -117.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -115.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -112.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -110.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -107.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -105.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -102.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -100.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -97.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -95.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -92.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -90.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -87.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -85.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -82.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -80.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -77.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -75.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -72.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -70.000000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -67.500000] Acquisition Number: 3, Image Position Patient: [-199.800003, -180.000000, -65.000000] Dimensiones del arreglo de segmentación: (300, 512, 512) Dimensiones del arreglo de imágenes CT: (194, 512, 512)
Confirmamos que hay más de una adquisición
Que hago aqui? Hay más de una adquisición pero solo debe haber una. En los obejtivos pone "Asegurar que solo haya una adquisicion".
In [ ]:
check_acquisition_uniformity(sorted_ct_files)
Multiple acquisition numbers found. Ensure all CT images are from a single acquisition.
Aplicamos windowing a las imagenes para mejorarlas¶
In [ ]:
# Windowing parameters
window_center, window_width = 40, 400
windowed_ct_images = np.array([apply_window(img, window_center, window_width) for img in ct_images])
In [ ]:
# Visualize different planes
visualize_planes(windowed_ct_images)
In [ ]:
# Visualize images and segmentations
for i in range(10): # Example for the first 10 slices
display_images(windowed_ct_images[i], segmentation_image[i])
print("Segmentation array loaded with shape:", segmentation_image.shape)
Segmentation array loaded with shape: (300, 512, 512)
ROI¶
In [ ]:
def apply_threshold(ct_slice, threshold=0.1):
"""Aplicar umbralización para separar el ROI del fondo."""
_, thresh = cv2.threshold(ct_slice, threshold, 255, cv2.THRESH_BINARY)
return thresh
'''
PREGUNTA: Tenemos que excluir aquellas zonasa que son demasiado pequeñas?
Se podría implementar un filtro por tamaño,
conservando solo las ROIs que tengan un número mínimo de píxeles,
lo que indica que son regiones significativas.
Ej. unique, counts = np.unique(labels_im, return_counts=True)
filtered_labels = np.zeros_like(labels_im)
for label, count in zip(unique, counts):
if count > min_pixel_count: # Filtrar ROIs pequeñas
filtered_labels[labels_im == label] = label
'''
def find_rois(thresh_image):
"""Encontrar y etiquetar ROIs en la imagen umbralizada."""
num_labels, labels_im = cv2.connectedComponents(thresh_image.astype(np.uint8))
return num_labels, labels_im
def visualize_rois(ct_slice, labels_im):
"""Visualizar las ROIs sobre la imagen original."""
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.imshow(ct_slice, cmap='gray')
plt.title('Original CT Slice')
plt.subplot(1, 2, 2)
plt.imshow(labels_im, cmap='nipy_spectral')
plt.title('Detected ROIs')
plt.show()
ct_slice = windowed_ct_images[0]
# Análisis de datos de la imagen
print("Estadísticas de la imagen:")
print(f"Valor mínimo: {np.min(ct_slice)}, Valor máximo: {np.max(ct_slice)}")
print(f"Media: {np.mean(ct_slice)}, Desviación estándar: {np.std(ct_slice)}")
# Procesamiento
threshold = 0.43
thresh = apply_threshold(ct_slice, threshold)
print(f"Porcentaje de píxeles sobre el umbral: {np.sum(thresh > 0) / thresh.size * 100:.2f}%")
num_labels, labels_im = find_rois(thresh)
print(f"Número de ROIs detectadas: {num_labels - 1}") # Excluyendo el fondo
# Visualización
visualize_rois(ct_slice, labels_im)
# Análisis de componentes conectados
'''
unique, counts = np.unique(labels_im, return_counts=True)
for roi_label, count in zip(unique, counts):
if roi_label != 0: # Excluyendo el fondo
print(f"ROI Label: {roi_label}, Pixel Count: {count}")
'''
Estadísticas de la imagen: Valor mínimo: 0.0, Valor máximo: 1.0 Media: 0.2165085220336914, Desviación estándar: 0.2696629915027511 Porcentaje de píxeles sobre el umbral: 25.47% Número de ROIs detectadas: 208
Out[ ]:
'\nunique, counts = np.unique(labels_im, return_counts=True)\nfor roi_label, count in zip(unique, counts):\n if roi_label != 0: # Excluyendo el fondo\n print(f"ROI Label: {roi_label}, Pixel Count: {count}")\n'
MIP¶
In [ ]:
def verify_alignment(ct_dcm, seg_dcm):
"""Verify if the DICOM files have matching slice locations."""
if hasattr(ct_dcm, 'SliceLocation') and hasattr(seg_dcm, 'SliceLocation'):
if abs(ct_dcm.SliceLocation - seg_dcm.SliceLocation) > 0.1: # Assuming a tolerance of 0.1
print(f"Mismatch: CT {ct_dcm.SliceLocation} vs. Seg {seg_dcm.SliceLocation}")
return False
print("CT and Segmentation are aligned.")
return True
print("One or both files do not contain SliceLocation data for comparison.")
return False
def visualize_ct_with_segmentation(ct_slice, seg_slice):
"""Visualize CT slice with segmentation overlay."""
fig, ax = plt.subplots(figsize=(10, 8))
ax.imshow(ct_slice, cmap='gray')
colored_cmap = create_transparent_cmap(plt.cm.viridis, 0.3)
ax.imshow(seg_slice, cmap=colored_cmap, interpolation='none')
ax.set_title('CT Slice with Segmentation Overlay')
plt.axis('off')
plt.show()
def create_transparent_cmap(cmap, transparency):
"""Create a colormap with transparency for overlay."""
colors = cmap(np.linspace(0, 1, 256))
colors[:, -1] = np.linspace(0, transparency, 256)
return matplotlib.colors.ListedColormap(colors)
def calculate_mip(ct_volume, axis):
"""Calculate the Maximum Intensity Projection (MIP) along a given axis."""
return np.max(ct_volume, axis=axis)
def calculate_and_display_mips(ct_volume):
"""Calculate MIPs for all three planes and display them with a lilac color theme."""
mip_axial = calculate_mip(ct_volume, axis=0)
mip_coronal = calculate_mip(ct_volume, axis=1)
mip_sagittal = calculate_mip(ct_volume, axis=2)
display_mip_images(mip_axial, mip_coronal, mip_sagittal)
def display_mip_images(mip_axial, mip_coronal, mip_sagittal):
"""Display MIP images for axial, coronal, and sagittal planes with customized color."""
plt.figure(figsize=(18, 6))
# Colormap personalizado que tiende al lila
cmap = plt.get_cmap('twilight')
# Definir los percentiles para ajustar los valores de visualización de cada MIP
vmin_a, vmax_a = np.percentile(mip_axial, [2, 98])
vmin_c, vmax_c = np.percentile(mip_coronal, [2, 98])
vmin_s, vmax_s = np.percentile(mip_sagittal, [2, 98])
# Axial
plt.subplot(1, 3, 1)
plt.imshow(mip_axial, cmap=cmap, aspect='auto', vmin=vmin_a, vmax=vmax_a)
plt.title('MIP Axial')
plt.axis('off')
# Coronal
plt.subplot(1, 3, 2)
plt.imshow(mip_coronal, cmap=cmap, aspect='auto', vmin=vmin_c, vmax=vmax_c)
plt.title('MIP Coronal')
plt.axis('off')
# Sagittal
plt.subplot(1, 3, 3)
plt.imshow(mip_sagittal, cmap=cmap, aspect='auto', vmin=vmin_s, vmax=vmax_s)
plt.title('MIP Sagittal')
plt.axis('off')
plt.tight_layout()
plt.show()
In [ ]:
# Example usage within your data processing pipeline:
'''for ct_dcm, seg_dcm in zip(ct_images, segmentation_image):
if verify_alignment(ct_dcm, seg_dcm):
visualize_ct_with_segmentation(ct_dcm.pixel_array, seg_dcm.pixel_array)'''
calculate_and_display_mips(ct_images)
3D¶
In [ ]:
def plot_3d(image, threshold=-300, rotation_angle=360, interval=50):
"""Plot 3D visualization of CT data."""
p = image.transpose(2, 1, 0)
p = p.clip(min=threshold)
fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection='3d')
pos = np.where(p > threshold)
scat = ax.scatter(pos[0], pos[1], pos[2], c=p[pos], cmap='viridis', s=1)
ax.view_init(30, 185)
ani = FuncAnimation(fig, lambda frame: ax.view_init(30, frame), frames=np.arange(0, rotation_angle, 2), interval=interval)
plt.colorbar(scat, ax=ax, orientation='vertical')
plt.show()
return ani
plot_3d(ct_images, threshold=100) # Adjust threshold as needed
Out[ ]:
<matplotlib.animation.FuncAnimation at 0x15f70622760>
Pregunta: Como quito la linea de atras? Esa columna que se ve cerca del eje y
Prueba: Superposición de las CT sobre la segmentada¶
In [ ]:
def visualize_ct_with_segmentation(ct_slice, seg_slice):
# Crear figura y ejes
fig, ax = plt.subplots(figsize=(10, 8))
ax.imshow(ct_slice, cmap='gray') # Visualizar la imagen de CT en escala de grises
# Configuración del colormap para la segmentación con transparencia
cmap = plt.cm.viridis
colors = cmap(np.linspace(0, 1, 256))
colors[:, -1] = np.linspace(0, 0.3, 256) # Configuración de la transparencia
colored_cmap = matplotlib.colors.ListedColormap(colors)
# Superponer la segmentación en la imagen de CT
seg_im = ax.imshow(seg_slice, cmap=colored_cmap, interpolation='none')
ax.set_title('CT Slice with Segmentation Overlay')
plt.axis('off')
plt.show()
# Iterar sobre cada slice en el conjunto de datos de CT y Segmentación
for i in range(5):
ct_slice = ct_images[i+5]
seg_slice = segmentation_image[i+5]
visualize_ct_with_segmentation(ct_slice, seg_slice)
Image Coregistration¶
Cargamos las Imagenes DICOM¶
In [ ]:
reference_image_path = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/Project/icbm_avg_152_t1_tal_nlin_symmetric_VI.dcm"
input_images_path = "C:/Users/Pedro/Documents/GitHub/Medial_Image_Pydicom/Project/RM_Brain_3D-SPGR"
In [ ]:
input_images = sort_dicom_files_by_position(load_dicom_files(input_images_path))
reference_image = pydicom.dcmread(reference_image_path)
In [ ]:
def multiply_quaternions(
q1: tuple[float, float, float, float],
q2: tuple[float, float, float, float]
) -> tuple[float, float, float, float]:
""" Multiply two quaternions, expressed as (1, i, j, k). """
return (
q1[0] * q2[0] - q1[1] * q2[1] - q1[2] * q2[2] - q1[3] * q2[3],
q1[0] * q2[1] + q1[1] * q2[0] + q1[2] * q2[3] - q1[3] * q2[2],
q1[0] * q2[2] - q1[1] * q2[3] + q1[2] * q2[0] + q1[3] * q2[1],
q1[0] * q2[3] + q1[1] * q2[2] - q1[2] * q2[1] + q1[3] * q2[0]
)
def conjugate_quaternion(
q: tuple[float, float, float, float]
) -> tuple[float, float, float, float]:
""" Multiply two quaternions, expressed as (1, i, j, k). """
return (
q[0], -q[1], -q[2], -q[3]
)
def mutual_information(hgram):
""" Calculate the mutual information between two images """
pxy = hgram / float(np.sum(hgram))
px = np.sum(pxy, axis=1)
py = np.sum(pxy, axis=0)
px_py = px[:, None] * py[None, :]
nzs = pxy > 0
return np.sum(pxy[nzs] * np.log(pxy[nzs] / px_py[nzs]))
def translation(
point: tuple[float, float, float],
translation_vector: tuple[float, float, float]
) -> tuple[float, float, float]:
""" Perform translation of `point` by `translation_vector`. """
x, y, z = point
v1, v2, v3 = translation_vector
return (x+v1, y+v2, z+v3)
def transform_point(point, t, q):
"""Apply translation and rotation to a point using quaternions."""
# Translate
point = np.array(point) + np.array(t)
# Rotate using quaternion
p = (0,) + tuple(point)
p_rotated = multiply_quaternions(multiply_quaternions(q, p), conjugate_quaternion(q))
return p_rotated[1:]
def axial_rotation(
point: tuple[float, float, float],
angle_in_rads: float,
axis_of_rotation: tuple[float, float, float]) -> tuple[float, float, float]:
""" Perform axial rotation of `point` around `axis_of_rotation` by `angle_in_rads`. """
x, y, z = point
v1, v2, v3 = axis_of_rotation
# Normalize axis of rotation to avoid restrictions on optimizer
v_norm = math.sqrt(sum([coord ** 2 for coord in [v1, v2, v3]]))
v1, v2, v3 = v1 / v_norm, v2 / v_norm, v3 / v_norm
# Quaternion associated to point.
p = (0, x, y, z)
# Quaternion associated to axial rotation.
cos, sin = math.cos(angle_in_rads / 2), math.sin(angle_in_rads / 2)
q = (cos, sin * v1, sin * v2, sin * v3)
# Quaternion associated to image point
q_star = conjugate_quaternion(q)
p_prime = multiply_quaternions(q, multiply_quaternions(p, q_star))
# Interpret as 3D point (i.e. drop first coordinate)
return p_prime[1], p_prime[2], p_prime[3]
def ssd(reference, transformed):
"""Compute Sum of Squared Differences between two images."""
return np.sum((reference - transformed) ** 2)
def transform_image(image, params):
"""Apply translation and rotation to the image based on params."""
# Unpack translation and rotation parameters
t1, t2, t3, angle, axis = params
translated = translation(image, (t1, t2, t3))
rotated = axial_rotation(translated, angle, axis)
return rotated
def objective_function(params, input_image, reference_image):
"""Objective function to minimize (SSD between reference and transformed input)."""
transformed_image = transform_image(input_image, params)
return ssd(reference_image, transformed_image)
def transform_image(image, params):
"""Apply translation and rotation to the image based on params."""
# Unpack translation and rotation parameters
t1, t2, t3, angle, axis = params
translated = translation(image, (t1, t2, t3))
rotated = axial_rotation(translated, angle, axis)
return rotated
def objective_function(params, input_image, reference_image):
"""Objective function to minimize (SSD between reference and transformed input)."""
transformed_image = transform_image(input_image, params)
return ssd(reference_image, transformed_image)
In [ ]:
# Initial transformation parameters: translation (tx, ty, tz) and rotation (angle, axis)
initial_params = [0, 0, 0, 0, (0, 0, 1)] # example initial guess
# Optimization to find the best parameters
result = minimize(
objective_function,
initial_params,
args=(input_images, reference_image),
method='BFGS'
)
# Apply the optimized transformation to the input image
optimized_params = result.x
transformed_image = transform_image(input_images, optimized_params)
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In[220], line 5 2 initial_params = [0, 0, 0, 0, (0, 0, 1)] # example initial guess 4 # Optimization to find the best parameters ----> 5 result = minimize( 6 objective_function, 7 initial_params, 8 args=(input_images, reference_image), 9 method='BFGS' 10 ) 12 # Apply the optimized transformation to the input image 13 optimized_params = result.x File c:\Users\Pedro\anaconda3\envs\AIV_P1\lib\site-packages\scipy\optimize\_minimize.py:530, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options) 51 def minimize(fun, x0, args=(), method=None, jac=None, hess=None, 52 hessp=None, bounds=None, constraints=(), tol=None, 53 callback=None, options=None): 54 """Minimization of scalar function of one or more variables. 55 56 Parameters (...) 528 529 """ --> 530 x0 = np.atleast_1d(np.asarray(x0)) 532 if x0.ndim != 1: 533 raise ValueError("'x0' must only have one dimension.") ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (5,) + inhomogeneous part.
In [ ]:
from scipy.optimize import least_squares
def registration_error(params, ref_points, input_points):
"""Calculate error between reference points and transformed input points."""
t = params[:3]
q = (params[3],) + tuple(params[4:7]) # Quaternion parameters
transformed_points = [transform_point(inp, t, q) for inp in input_points]
return np.linalg.norm(np.array(transformed_points) - np.array(ref_points), axis=1)
# Puntos de referencia random
ref_landmarks = np.array([[10, 20, 30], [40, 50, 60], [70, 80, 90]])
input_landmarks = np.array([[12, 22, 32], [42, 52, 62], [72, 82, 92]])
# Initial parameters: [tx, ty, tz, qw, qx, qy, qz]
initial_params = [0, 0, 0, 1, 0, 0, 0] # Identity quaternion and no translation
result = least_squares(registration_error, initial_params, args=(ref_landmarks, input_landmarks))
print("Optimal parameters:", result.x)
Optimal parameters: [-0.01442782 -0.01436382 -0.01429953 0.98456473 0.00100211 -0.00199494 0.000993 ]